• Felicia Magpantay, Matthieu Domenech de Cellès, Carles Breto, Jennie Lavine, Daihai He
• Ed Ionides, Pej Rohani, Mercedes Pascual
• Alex Smith
• US National Institutes of Health
• US NIH/DHS Research and Policy for Infectious Disease Dynamics Program
• US National Center for Ecological Analysis and Synthesis & the “Inference for Mechanistic Models” Working Group

Let \(b\) be biology and \(m\), mathematics.
Then \(b\) is more interesting than \(m\).

Proof

\(m\) is an activity of \(M\), the set of mathematicians. \(M\subset B\), where \(B\) is the set of biological organisms. Therefore, \(b\) is at least as interesting as \(m\). Moreover, there is at least one element of \(b\) that has not yet been mathematically studied. Q.E.D.

In seeking to understand biological phenomena, models are tools for looking at the data.

They focus attention on the discrepancies between our hypotheses and reality.

How to quantify these discrepancies?

The key principle of statistics is that we must model the error.

• In seeking to infer causal mechanism in intact biological systems, dynamics are especially informative.
• Hypotheses about causal mechanisms can often readily be formalized as dynamical systems models.
• The problems of inference
  • How can we estimate unknown parameters?
  • How can we decide among competing models?
  • How can we extract maximum information from time series data?
  • When do we know when we need a better model vs better data?

The key principle is that we must model the error.

• Measurement error
• Process noise
  • Environmental trends and fluctuations
  • Unmodeled heterogeneities
  • Model misspecification
  • Unmodeled variables

1. Markov models
2. Examples: cholera and pertussis
3. Partially observed Markov processes
4. Plug and play inference methods
5. Plug and play methods in practice
6. Conclusions and recommendations

\[ \begin{aligned} dN_{XY}&=\lambda\,X\,dt+\sigma_{XY}(X)\,dW_1(t)\\ dN_{X\emptyset}&=\mu\,X\,dt+\sigma_{X\emptyset}(X)\,dW_2(t)\\ \end{aligned} \]

\[dW_1,dW_2\;\stackrel{i.i.d}{\sim}\;\mathrm{Normal}(0,\sqrt{dt})\]

\[ \begin{aligned} dX&=B\,dt-dN_{XY}-dN_{X\emptyset}\\ &=B\,dt-(\lambda+\mu)\,X\,dt\\ &\qquad-\sigma_{XY}(X)\,dW_1(t)-\sigma_{X\emptyset}(X)\,dW_2(t)\\ dY&=dN_{XY}\\ &=\lambda\,X\,dt+\sigma_{XY}(X)\,dW_1(t)\\ \end{aligned} \]

1. Markov models
2. Examples: cholera and pertussis
3. Partially observed Markov processes
4. Plug and play inference methods
5. Plug and play methods in practice
6. Conclusions and recommendations

1. Relative importance of human-human vs environmental transmission
2. Duration of infection-induced immunity
3. Roles of seasonal and multiennial climate drivers
   • complex seasonality
   • multiennial cycles
4. Importance of bacteriophage in environment

• nonlinear
• stochastic
• nonstationary
• extrinsically driven

• How well do each of the models explain the data?
• What insights can we gain into the system’s dynamic self-regulation?

Note that the structure of the models corresponds to the questions of scientific interest, without regard to issues of mathematical tractability

1. Relative importance of human-human vs environmental transmission
2. Duration of infection-induced immunity
3. Roles of seasonal and multiennial climate drivers
   • complex seasonality
   • multiennial cycles
4. Importance of bacteriophage in environment

(Lavine, King, and Bjørnstad 2011)

(Domenech de Cellès, Magpantay, King, and Rohani 2016)

Why the resurgence?

• Changes in vaccine efficacy
• Vaccine-driven pathogen evolution
• Increased circulation of congeneric pathogens
• Asymptomatic transmission
• Loss of natural immune boosting
• Increased surveillance sensitivity

Key unknowns:

• durability of vaccine-induced immunity
• relative efficacy of natural- and vaccine-derived immunity

\[\lambda(t) = \frac{\beta(t)\,(I_1+\theta\,I_2)+\bar{\beta}\,\iota}{N}\]

Post-vaccination infections are observed at a reduced rate, \(\eta\).

(Magpantay, Domenech de Cellès, Rohani, and King 2016)

• Can we estimate the nature and durability of vaccine-induced protection?
• Which models explain the data adequately?
• Questions of parameter estimation and model selection

(Magpantay et al. 2016)

(Domenech de Cellès, Magpantay, King, and Rohani 2018)

(Domenech de Cellès et al. 2018)

(Domenech de Cellès et al. 2018)

(Domenech de Cellès et al. 2018)

(Domenech de Cellès et al. 2018)

(Domenech de Cellès et al. 2018)

(Domenech de Cellès et al. 2018)

(Domenech de Cellès et al. 2018)

(Domenech de Cellès et al. 2018)

(Domenech de Cellès et al. 2018)

(Domenech de Cellès et al. 2018; Rohani, Zhong, and King 2010)

(Fumanelli, Ajelli, Manfredi, Vespignani, and Merler 2012; Mossong et al. 2008)

• Can we estimate the nature and durability of vaccine-induced protection?
• Which models explain the data adequately?
• Questions of parameter estimation and model selection

(Domenech de Cellès et al. 2018)

1. Markov models
2. Examples: cholera and pertussis
3. Partially observed Markov processes
4. Plug and play inference methods
5. Plug and play methods in practice
6. Conclusions and recommendations

(King, Nguyen, and Ionides 2016)

(King et al. 2016)

(King et al. 2016)

(King et al. 2016)

(King et al. 2016)

• observation times: \(t_1 < t_2 < \dots < t_N\)
• data: \(y_{t_n}^*\)
• unobserved Markov state process: \(X_t\)
• measurement process: \(Y_t\)
• process model (Markov state transition density): \(f_n(X_{t_n}|X_{t_{n-1}},\theta)\)
• measurement model: \(g_n(Y_{t_n}|X_{t_{n-1}},\theta)\)
• initial state distribution \(f_0(X_{t_0},\theta)\)

(King et al. 2016)

The full joint density is:

\[f_{X,Y}(x,y;\theta) = f_0(x_0;\theta)\,\prod_{n=1}^N\!f_{n}(x_n|x_{n-1};\theta)\,g_{n}(y_n|x_n;\theta).\]

The likelihood function is the marginal density for \(Y\), evaluated at the data:

\[ \begin{split} \mathcal{L}(\theta)&=f_{Y}(y^*_1,\dots,y^*_n;\theta)\\ &=\int f_{X,Y}(x_0,\dots,x_N,y^*_1,\dots,y^*N;\theta)\, dx_1\cdots dx_n. \end{split} \]

• computing the likelihood requires integrating over \(X_t\)
  for all \(t\)
• traditional solution—be very clever:
  • simplify the model to make it statistically tractable
  • devise and tune an inference algorithm (e.g., MCMC sampler) specifically to the model
• but: are we testing just our hypothesis or the statistical method also?
• it is perilous to invest too much time in one model
• the Pygmalion problem
• an antidote: the method of multiple working hypotheses

Definition: an algorithm has the plug-and-play property if it has no need to evaluate the process model density \(f(X_t|X_{t-1})\).

Plug-and-play methods access the process model only via simulation.

This puts minimal restrictions on the form of the models that can be entertained.

(He, Ionides, and King 2010; King et al. 2016)

1. Markov models
2. Examples: cholera and pertussis
3. Partially observed Markov processes
4. Plug and play inference methods
5. Plug and play methods in practice
6. Conclusions and recommendations

• Feature-based methods
  • Approximate Bayesian computation (ABC)
  • Nonlinear forecasting (NLF)
  • Attractor-reconstruction based methods
  • Probe matching
  • Synthetic likelihood
• Full-information methods (likelihood based)
  • Sequential Monte Carlo (the particle filter)
  • Particle Markov chain Monte Carlo
  • Iterated filtering

Available in the pomp software package (King et al. 2016): [https://kingaa.github.io/pomp/]

Given parameters \(\theta\in\Theta\) and data \(y^*\in\mathcal{Y}\)

• Select probes (summary statistics), \(\mathbb{S}:\mathcal{Y}\to\mathbb{R}^d\)
• Compute observed probes, \(s^* = \mathbb{S}(y^*)\)
• Simulate \(J\) datasets \(y_j\) and compute simulated probes \(s_j=\mathbb{S}(y_j)\)
• Using sample mean. \(\mu = \langle s_j \rangle\), and variance, \(\Sigma = \langle (s_j-\mu)\,(s_j-\mu)^T \rangle\), compute the synthetic log likelihood \[\ell_{\mathbb{S}}(\theta) = -\tfrac{1}{2}\,(s^*-\mu)^{T}{\Sigma}^{-1}(s^*-\mu)-\tfrac{1}{2}\,\log((2\pi)^d|\Sigma|)\]
• E.g., maximize \(\ell_{\mathbb{S}}\) over \(\Theta\)

(Fasiolo, Pya, and Wood 2014; Wood 2010)

• Given \(\theta\) and \(y^*\), initialize particle swarm by drawing \(X^F_{t_0,j}\sim f_0(\cdot|\theta)\), \(j=1,\dots,J\)
• For \(n = 1,\dots,N\),
  • Simulate for prediction: \(X^P_{t_{n},j}\sim f_k(\cdot|X^F_{t_{n-1},j},\theta)\)
  • Evaluate weights: \(w_{n,j} = g(y^*_{t_{n}}|X^P_{t_{n},j},\theta)\)
  • Perform weighted resample of \(\{X^P_{t_n,j}\}_{j=1}^J\) to obtain \(\{X^F_{t_n,j}\}_{j=1}^J\)
  • Estimate conditional log likelihood: \(\ell_{n}(\theta) = \log\mathrm{Prob}[y^*_{t_n}|y^*_{t_1},\dots,y^*_{t_{n-1}}] = \log\tfrac{1}{J}\sum_j\!w_{n,j}\)
• Unbiased estimate of the log likelihood: \(\ell(\theta) = \sum_n\,\ell_n(\theta)\)

(Arulampalam, Maskell, Gordon, and Clapp 2002)

• Given \(\theta_0\), compute \(\ell(\theta_0)\) using the PF
• For \(m>0\), draw \(\theta^P_m\sim q(\cdot|\theta_{m-1})\)
• Compute \(\ell(\theta^P_m)\) using the PF
• Generate \(U\sim\mathrm{Uniform}(0,1)\)
• Set \(\theta_m=\begin{cases}\theta^P_m, &\text{if } U<\frac{f_\Theta(\theta^P_m)\exp(\ell(\theta^P_m))}{f_\Theta(\theta_{m-1})\exp(\ell(\theta_{m-1}))},\\ \theta_{m-1},&\text{otherwise.}\end{cases}\)
• Repeat until the Markov chain converges

Because the particle filter gives as unbiased, though noisy, estimate of the likelihood, this MCMC converges to the correct target distribution.

(Andrieu, Doucet, and Holenstein 2010)

• main problem with particle filter: particle depletion
• basic idea of iterated filtering:
  • replace fixed parameters by parameters that execute a random walk, e.g., \[\theta \mapsto \theta_n \sim \mathrm{Normal}(\theta_{n-1},\sigma)\]
  • this combats particle depletion by introducing variability
  • also permits local exploration of the likelihood surface
  • take \(\sigma\to 0\) iteratively
• analogies with simulated annealing

(Ionides, Bretó, and King 2006; Ionides, Nguyen, Atchadé, Stoev, and King 2015; King et al. 2016)

• Choose starting point \(\theta_0\) and random-walk perturbations \(\sigma_{m}\), \(m=0,1,2,\dots\).
• Initialize particle swarm by drawing \(\theta^F_{j,0} \sim \mathrm{Normal}(\theta_0,\sigma_0)\)
• For \(m>0\), run PF to obtain \(\{\theta^F_{j,m}\}_{j=1}^J\)
• Repeat unto convergence
• The \(\theta^F\) converge to the maximum likelihood estimate of \(\theta\)

(Ionides et al. 2015)

1. Markov models
2. Examples: cholera and pertussis
3. Partially observed Markov processes
4. Plug and play inference methods
5. Plug and play methods in practice
6. Conclusions and recommendations

• The best model agrees with already published results
• The best model is the one that I can estimate before my thesis is due
• The best model is the one that best agrees with my preconceptions
• The best model is the one that makes my supervisor the happiest
• An identifiable model is better than a non-identifiable model
• The best model is the one with the narrowest confidence intervals

• How well does the model fit the data?
  • How does the likelihood compare to benchmark values?
• Fitting best \(\ne\) fitting well
• Does the model fit well?
  • What aspects of the data are explained?
  • What aspects are ascribed to noise?
• Model adequacy
  • Are the data consistent with the model’s assumptions?
  • If the model were correct, how plausible are the data?

• Predictive power
  • Can the model predict new data?
  • What predictions to parameter estimates make about the results of different kinds of studies?
• Fertility
  • Does the model make testable predictions?

(King, Ionides, Pascual, and Bouma 2008)

(King et al. 2008)

(King et al. 2008)

Is it necessary that all infected individuals be equally infectious?
(King et al. 2008)

case fatality: \(0.34\pm 0.2\)
duration of immunity \(1.5 \pm 0.7~\text{yr}\)
(King et al. 2008)

\[\lambda(t) = \frac{\beta(t)\,(I_1+\theta\,I_2)+\bar{\beta}\,\iota}{N}\]

Post-vaccination infections are observed at a reduced rate, \(\eta\).

(Magpantay et al. 2016)

Interpretation: in vaccinated hosts, infections are mild to asymptomatic, yet equally infectious

(Magpantay et al. 2016)

Flat profiles indicate lack of information in the data relative to the question

(Magpantay et al. 2016)

• relative transmissibility 0.99 [0.40–1.00]
• relative observability 0.39 [0.19–1.00]
• high detection rates in adolescents and in adults: 24% [10–66%]
• despite effectiveness, eradication via routine immunization alone is not possible
  (Domenech de Cellès et al. 2018)

1. Markov models
2. Examples: cholera and pertussis
3. Partially observed Markov processes
4. Plug and play inference methods
5. Plug and play methods in practice
6. Conclusions and recommendations

• SL throws away information, but gains robustness
• Many choices are possible for SL probes, with consequences for inference
• IF converges much faster than PMCMC
• Frequentist vs Bayesian inference styles
• Smoothed state trajectories can be computed by PMCMC
• Probes can be used for post hoc checks of model adequacy

• Use IF, profile likelihood, and profile trace plots to identify heights in the likelihood surface.
• Use PMCMC to obtain posterior densities and smoothed state trajectories if desired.
• Use independent parameter estimates to diagnose model misspecification.
• Use probes to ascertain model adequacy and diagnose model misspecification.

• It is essential to model the error, including both process noise and measurement error.
• Plug-and-play inference methods facilitate scientific investigation.
• Effective, statistically efficient plug-and-play inference methods are available.
• There is an intense need for further methodological development of such methods to improve computational efficiency and accommodate new data types.